Numerical Linked-Cluster Algorithms. II. t-J models on the square lattice 
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We discuss the application of a recently introduced numerical linked-cluster (NLC) algorithm to 
strongly correlated itinerant models. In particular, we present a study of thermodynamic observ- 
ables: chemical potential, entropy, specific heat, and uniform susceptibility for the t-J model on the 
square lattice, with J/t = 0.5 and 0.3. Our NLC results are compared with those obtained from 
high-temperature expansions (HTE) and the finite-temperature Lanczos method (FTLM) . We show 
that there is a sizeable window in temperature where NLC results converge without extrapolations 
whereas HTE diverges. Upon extrapolations, the overall agreement between NLC, HTE, and FTLM 
is excellent in some cases down to 0.25t. At intermediate temperatures NLC results are better con- 
trolled than other methods, making it easier to judge the convergence and numerical accuracy of 
the method. 

PACS numbers: 05.50.+q,05.70.-a,75.10.Lp,05.10.-a 
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I. INTRODUCTION 

In a recent paper [l| we introduced a linked-cluster 
algorithm, which we called the numerical linked-cluster 
(NLC) algorithm, that allows one to obtain temperature- 
dependent properties of quantum lattice models in the 
thermodynamic limit from the exact diagonalization of 
small clusters. A detailed exposition of NLC and its ap- 
plication to quantum spin models on square, triangular, 
and kagome lattices has been presented in Ref. @ . There 
we have shown that for many spin models and thermody- 
namic quantities NLC results can be substantially more 
accurate than HTE and exact diagonalization (ED). 

In this paper we discuss how to use NLC to calcu- 
late properties of strongly correlated itinerant models. 
In particular, we study the thermodynamics of the pla- 
nar t-J model on the square lattice. This model was 
introduced by Anderson and others @, 0] as a means 
to understanding the microscopic mechanism for high- 
temperature superconductivity @, Q. It is one of the 
simplest models that allows one to study the interplay 
between itinerancy of electrons and their spin fluctua- 
tions, possibly leading to superconductivity and many 
other exotic quantum phases Q- 

In spite of its simplicity, understanding finite- 
temperature thermodynamic properties of the t-J model 
has proven to be a very challenging task Q|. Quantum 
Monte Carlo simulations suffer from severe sign prob- 
lems, which become a major difficulty at low tempera- 
tures. The two general approaches that have been com- 
monly used to study this model arc ED and HTE. ED 
studies in which one fully diagonalizes the t-J Hamilto- 
nian are difficult since they can only be done for very 
small systems [1, as a consequence of which fi- 

nite size effects are very large. A more efficient ap- 
proach to this problem is the finite-temperature Lanczos 
method (FTLM), which has been developed by Jaklic 



and Prelovsek (JP) [ri| . Within this approach the full 
thermodynamic trace is reduced by randomly sampling 
the eigenstates of the Hamiltonian. This allows one to 
study larger systems sizes in an unbiased way, but still 
finite size effects become relevant as the temperature is 
lowered. 

In order to obtain results in the thermodynamic 
limit one can use high-temperature expansions (HTE) 
[iH . [H, [2 Within this method the properties of 

the system are expanded in powers of the inverse tem- 
perature (3 0, [ijj • These expansions, carried out to 
order /3 N , where N is of order 10, provide accurate nu- 
merical results within the radius of convergence of the 
series. The temperature scale where HTE converges can 
be rather high for t-J models. It is typically set by t 
when t > J. Beyond the region of convergence, series 
extrapolation methods allow one to calculate ther- 
modynamic properties, but their reliability remains un- 
certain. Because the region of convergence is small in in- 
verse temperature, extrapolations even to temperatures 
of order J are more sensitive to the choice of the extrap- 
olation method and variables and hence less reliable than 
for purely spin models. 

Here, we study thermodynamic properties of the t-J 
model using NLC. The basic idea is to show with exam- 
ples the advantages and disadvantages of NLC as com- 
pared with FTLM and HTE, and the region of tempera- 
tures that within NLC can be accessed for the different 
observables of interest. The exposition is organized as 
follows. In Sec. [IT] we introduce the t-J Hamiltonian and 
discuss the basic ideas of the NLC calculation for itiner- 
ant models. The calculations must be done in the grand 
canonical ensemble and the change from chemical poten- 
tial to density must be done numerically. We will show 
that this works quite well. In the remaining sections we 
present results for the entropy (Scc. lIIl| ). the uniform sus- 
ceptibility (Sec. IIV|) . and the specific heat (Sec.|V|. The 
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conclusions are presented in Sec. I VII 



II. GRAND CANONICAL NLC 



The t-J Hamiltonian can be written as 



n 



(Si ■ Sj - -mrij 



where c\ s and c is are the creation and annihilation op- 
erators for an electron with spin s =1,1 on a site i, 
rii = c\ s c is is the density operator, P is a projection 
operator to ensure no hopping produces doubly occupied 
sites, i.e., we assume that the local Coulomb repulsion is 
very large such that two electrons (with antiparallel spin) 
cannot be on the same lattice site, and 



S« — 2 E/ C is (J ss'C i 



(2) 



is the local spin operator (a are the Pauli matrices). The 
sums in Eq. ([T]) run over nearest-neighbor sites. 

As in the case of spin systems [H the fundamen- 
tal basis of our numerical linked cluster expansion, for 
some extensive property P of an infinite lattice C, is the 
relation 111, [T3l 



P(C)/N = J^L(c)xW P (c), 



In Eq. ([5]), /i denotes the chemical potential, T is the 
temperature of the system (we have set the Boltzmann 
constant ks to be unity), and Z c is the partition function 
in each cluster c 



Z c , = Tr 



exp 



IT 



(0) 



Notice that in Eq. ([5]) and Eq. ([6]), Tr denotes the grand 
canonical trace, i.e., we fully diagonalize the t-J Hamil- 
tonian for all possible fillings in each cluster. This, to- 
(1) gcther with the fact that the Hilbcrt space is larger (three 
states per site) for the t-J model, when compared with 
1/2-spin models (with only two states per site) studied 
in Refs. [HQ, reduces the size of the largest clusters we 
can consider here. All results presented in this paper cor- 
respond to the site expansion of the square lattice (see 
Ref. for details) with all clusters up to 10 sites. At 
zero doping, i.e., the Heiscnbcrg model, we present NLC 
results considering all clusters with up to 13 sites [H, 0]- 
In addition, throughout this work we examine two values 
of the parameter J/t: J/t = 0.5 as in the HTE studies of 
RefJIl, and J/t = 0.3 as in the FTLM studies of Refs. 
fiEITto both of which we compare some of our results. 

For itinerant models, the first quantity one needs to 
evaluate within NLC is the density. This is because one 
is, in general, interested in the behavior of thermody- 
namic observables as a function of the temperature at 
a fixed density or as a function of the density (or hole 
doping) at a fixed temperature. In Fig. [1] we plot the 
hole density as a function of the temperature for differ- 
ent chemical potentials. The results of the bare NLC 
(o\ sums exhibit a clear feature. They converge to lower 
temperatures as the density approaches 1. 



where the left hand side is the value of the property P 
per lattice site in the thermodynamic limit. On the right 
hand side L(c) is the so-called lattice constant that is the 
number of embeddings of the cluster c, per lattice site, 
in the lattice C. Wp(c) is the weight of the cluster c for 
the property P. The latter is defined recursively by the 
principle of inclusion and exclusion [l(| 



W P (c)=P(c)-Y^Wp(i 



(4) 



where P(c) is the property P calculated for the finite 
cluster c. The sum on s runs over all subclusters of c. 
The basic idea of NLC is to calculate P and Wp at any 
temperature by means of an exact diagonalization of each 
cluster c. 

For itinerant models, it is desirable to control the den- 
sity in the thermodynamic limit. Within NLC this is 
achieved by working in the grand canonical ensemble, 
i.e., by introducing a chemical potential. Hence, in this 
case 



P(c) = ^Tr \ P exp 



(5) 




FIG. 1: (Color online) Hole density as a function of temper- 
ature for different chemical potentials (from top to bottom 
^ = 0, 0.2, . . . , 2.4), and J/t = 0.5. For each value of the 
chemical potential we have plotted results of the bare NLC 
sums up to 9 and 10 sites (thin continuous lines) and results of 
Wynn extrapolations after three and four cycles of improve- 
ment (see Ref. @] for details). 

In order to access lower temperatures we use the se- 
quence extrapolation techniques detailed in Ref. Q . Re- 
sults of Wynn extrapolations for the density after three 
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and four cycles of improvement are also depicted in Fig. 
[TJ For the lowest chemical potential (fi = 0) one can 
see that series extrapolations allow one to extend the re- 
gion of convergence of the density from T/t > 0.75 to 
T/t < 0.2. 

Once we have performed the extrapolations for the 
density, we numerically invert the dependence of any ob- 
servable from (//, T) to (n, T). As a first example we show 
in Fig.[2]the dependence of the chemical potential on the 
temperature when the density is held fixed. As can be 
seen, it is possible to follow the chemical potential curves 
at constant density to rather low temperatures T/t ~ 0.1. 
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FIG. 2: (Color online) Chemical potential as a function 
of temperature for different densities (from bottom to top 
n = 0.7, 0.72, . . . , 0.98), and J/t = 0.5. For each value of the 
density we have plotted results of Wynn extrapolations after 
three and four cycles of improvement. 

An analysis of the results depicted in Fig. [2] for the 
lowest temperatures suggests that the dependence of /x 
on the temperature is roughly linear 

M (T) = n(T = 0) + A T. (7) 

This is in agreement with the results obtained by JP @ 
using FTLM, and departs from the T 2 dependence ex- 
pected of a Fermi liquid. In addition, there is a change 
in the slope A as the density is shifted. For J/t = 0.5 
we find that the slope A changes sign when n! ~ 0.81. 
This can be better seen in Fig. [3j where we plot the hole 
density vs the chemical potential at fixed temperatures. 
There, one can see that for 1 — n' ~ 0.19 the different 
curves cross each other. For 1—n' < 0.19 our NLC results 
are well converged at all temperatures T/t > 0.15. Small 
departures between Wynn extrapolations after three and 
four cycles of improvement are only apparent for the low- 
est temperatures and highest hole concentrations. 

We have also performed a similar study for J/t = 0.3. 
In this case we have obtained, in close agreement with JP 
Q, that the change of slope in A occurs for n' ~ 0.86. 




FIG. 3: (Color online) Hole density vs chemical potential for 
different temperatures (below the crossing point from bottom 
to top T/t = 0.15, 0.20, ...,0.5), and J/t = 0.5. For each 
value of the temperature we have plotted results of Wynn 
extrapolations after three and four cycles of improvement. 

III. ENTROPY 

Once one has numerically inverted the dependence 
(/i, T) to (n,T) it is possible to study the behavior of 
observables of interest as a function of the temperature 
at a fixed density. In this section we consider the entropy 
per lattice site 

S^(AnZ + ^-^ ), (8) 

where the term proportional to the chemical potential 
on the right hand side, not present in our calculations 
for spin systems, is needed when dealing with the grand 
canonical ensemble. ./V is the number of lattice sites. 

In Fig. |U we plot the entropy vs temperature for dif- 
ferent values of the density after three and four cycles 
of improvements using Wynn's algorithm. As with the 
chemical potential, the NLC convergence for S improves 
as the density approaches 1 . For all the densities plotted 
in Fig. 5] we get good convergence for T/t > 0.25. 

Now we discuss the dependence of the entropy on den- 
sity when the temperature is held fixed. This is shown in 
Fig. [5j The entropy exhibits a very broad maximum that 
slowly shifts towards lower hole densities as the temper- 
ature is increased. These results are in qualitative agree- 
ment with experimental measurements that have found 
the entropy to be maximum around 1 — n = 0.22 at a 
temperature T/t = 0.07 (lower than the ones we have 
calculated here) [l9| . 

To conclude this section on the entropy we compare in 
Fig. [HI our NLC results with those obtained by JP Q us- 
ing FTLM in clusters with 20 lattice sites and J/t = 0.3. 
At zero doping (Heisenberg model), the agreement be- 
tween NLC (with clusters up to 13 sites Q) and FTLM 
is remarkable down to T/t ~ 0.25. Below that temper- 
ature, finite size effects start to be apparent as FTLM 
results marginally depart from the ones obtained using 
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FIG. 4: (Color online) Entropy as a function of temperature 
for different densities and J/t = 0.5. For each value of the 
density we have plotted results of Wynn extrapolations after 
three and four cycles of improvement. 



0.5 




FIG. 5: (Color online) Entropy as a function of the den- 
sity for different temperatures (from bottom to top T/t — 
0.25, 0.20, . . . , 0.6) and J/t = 0.5. For each value of the tem- 
perature we have plotted results of Wynn extrapolations after 
three and four cycles of improvement. 



NLC (see also discussion in Ref. [l| on the entropy of the 
Heisenberg model on the square lattice). 

Within FTLM, finite size effects for the entropy, at 
a given temperature, increase when holes are added to 
the antiferromagnct. This can be understood within the 
analysis presented in Ref. Q. There it was shown that 
the temperature at which finite size effects start to be- 
come relevant increases with doping 1 — n > 0.15, for the 
system sizes considered in that work. In particular, for 
n = 0.7 (the worst case for FTLM) our NLC results for 
S are slightly different from the ones in Ref. Q for all 
temperatures in Fig. [6l 



IV. UNIFORM SUSCEPTIBILITY 

We study in this section another thermodynamic quan- 
tity of much experimental interest, the uniform suscepti- 



1.0 - 



to 



0.5 



0.0 



n=0.7 - 



n=1.0 




FTLM (JP) 
NLC Wynn 4 



0.0 



0.5 



T/t 



1.0 



1.5 



FIG. 6: (Color online) Entropy as a function of temperature 
for different densities (from top to bottom at T/t = 1.5 NLC 
results correspond to n — 0.7, 0.8, 0.9, and 1.0) and J/t = 0.3. 
Our Wynn extrapolations after four cycles of improvement are 
compared with FTLM results by JP @] . 
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In Fig. [7]we show results for the uniform susceptibility 
as a function of temperature for different densities. For 
the Heisenberg model, the results of our Wynn extrapola- 
tions for x are well converged down to T / J ~ 0.3. They 
allow one to resolve the peak in x that occurs around 
T/J ~ 1, and signals the onset of short range antifcrro- 
magnctic order. With doping the peak shifts to lower 
temperatures, i.e., the doping slows down the growth 
of antiferromagnetic correlations, and eventually disap- 
pears. Only for n = 0.9 and 0.95 do our NLC results 
exhibit signals of a peak in x- 
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FIG. 7: (Color online) Uniform susceptibility as a function 
of the temperature for different densities (from bottom to 
top for the highest temperatures n = 0.8, 0.85, . . . , 1.0) and 
J/t = 0.5. For each value of the density we have plotted 
results of Wynn extrapolations after three and four cycles of 
improvement. 
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In what follows we compare NLC results for the uni- 
form susceptibility with the ones obtained using HTE 
flU ). We first contrast, in Fig. the results of the NLC 
bare sums ([3]) for the bond and site expansions with those 
obtained for the HTE bare sums up to order 10. Fig- 
ure [U shows that the direct NLC sums converge down to 
T/t <~ 0.75 (with the site expansion being slightly bet- 
ter than the bond expansion) , while the HTE results are 
only well converged to T/t ~ 1.25. The existence of this 
region of temperatures where NLC converges while HTE 
does not was our main motivation for developing NLC, 
and makes NLC a more controlled technique at interme- 
diate temperatures. 
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FIG. 8: (Color online) Uniform susceptibility results for the 
bare NLC sums of the bond and site expansions compared 
with those obtained with HTE up to order 10. The density 
in the system is n = 0.85. 



As discussed before, an important feature of NLC and 
HTE is that both approaches allow for systematic ex- 
trapolations that accelerate the convergence of NLC and 
enable going beyond the radius of convergence of HTE. In 
Fig. [5]wc compare Padc extrapolations for x with Wynn 
extrapolations for the NLC site expansion. The agree- 
ment between these two approaches is remarkable down 
to T/t r~j 0.5, which show that indeed extrapolations can 
work quite well for both techniques, if the temperature is 
not too low. For this quantity, however, different Wynn 
extrapolations converge to lower temperatures than the 
Pade approximants. We find this to be remarkable for 
itinerant models because the analytic structure of HTE 
allows for an analytic inversion of the grand-canonical 
(/i, T) dependence into the more experimentally relevant 
one (n, T). Within NLC we have to perform a numerical 
extrapolation to access lower temperatures and then do 
a numerical inversion (/i, T)— >(n, T). Still, the fact that 
NLC contains exact information from the finite clusters 
at all temperatures, while in HTE only a power series 
expansion in inverse temperature is kept, means that the 
former, in a real sense, requires less extrapolation. 
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FIG. 9: (Color online) Pade approximants for the uniform 
susceptibility [l3|] are compared with Wynn extrapolations of 
the NLC site expansion. In the legend, Pade m, Pade u, and 
Pade 1, indicate the mean estimated value, upper, and lower 
limits, respectively, obtained by different Pade approximants. 
The density in the system is n = 0.85. 



V. SPECIFIC HEAT 

In this section wc analyze the NLC results for the spe- 
cific heat and compare them with those obtained within 
FTLM The specific heat is defined as 



-as 

' dT 



(10) 
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However, we avoid the numerical differentiation of the 
entropy with respect to the temperature by evaluating 
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which substantially reduces the numerical errors as the 
NLC sums for the number of particles (N e = (X^ n »)) 
and the energy ((H)) converge (and are better behaved) 
to lower temperatures than the ones of the entropy. Still, 
we have to perform numerical derivatives in addition to 
calculating the fluctuations of the energy ((H 2 ) — (H) 2 ). 
Hence, of the thermodynamic obscrvablcs analyzed so far 
C v is the most difficult for NLC evaluation. 

In Fig.[Tn]we plot the specific heat as a function of the 
temperature for different densities. For C v we only obtain 
well converged results below T/t = 0.5 for the Heisen- 
berg model and the t-J model at low hole concentration. 
Figure fTDl shows that the maximum in the specific heat, 
attributed to the thermal activation of the spin degrees of 
freedom, becomes strongly suppressed with doping. This 
can be better seen in Fig. [11] where we have plotted the 
specific heat as a function of the hole density for fixed 
values of the temperature. For the lowest temperatures 
in that figure, we only can follow the C v curves for doping 
below 0.15. They clearly show that C v decreases as the 
doping is increased, and exhibits a minimum that moves 
towards higher doping concentrations as the temperature 
is lowered. 
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FIG. 10: (Color online) Specific heat as a function of the 
temperature for different densities (from top to bottom n = 
0.7, 0.75, . . . , 1.0) and J/t = 0.5. For each value of the density 
we have plotted results of Wynn extrapolations after three 
and four cycles of improvement. 




FIG. 11: (Color online) Specific heat as a function of the 
hole density for different temperatures (from top to bottom 
n = 0.5, 0.55, . . . , 1.0) and J/t = 0.5. For each value of the 
temperature we have plotted results of Wynn extrapolations 
after three and four cycles of improvement. 



In Ref. [2| we have argued that the specific heat of 
the Hcisenbcrg model suffers from large finite size effects 
at relatively high temperatures. (We presented results 
for the full diagonalization of a 4 x 4 lattice.) In Fig. 
[T2]we compare our NLC results for the t-J model with 
the ones obtained by JP for clusters with 20 lattice sites 
and J/t = 0.3. At high temperatures T/t > 0.75 our 
NLC results are in very good agreement with the ones 
of JP. NLC results for C v , however, do not converge be- 
low T = 0.5 when the doping in the system is large. As 
the density approaches 1, and the convergence of NLC 
moves to lower temperatures, our linked-cluster results 
start to depart from the ones obtained with FTLM. In 
general, we find our calculation of the specific heat in the 
thermodynamic limit to be below the FTLM results for 
finite clusters. This deviation becomes particularly large 
for the Heisenberg model. The FTLM peak for C v when 
N = 20 occurs at larger temperatures than the one ob- 



tained within NLC. Interestingly, increasing the system 
size up to N — 26 does not help since the peak remains 
at the same position while becoming higher Q, instead 
of moving to lower temperatures and becoming smaller 
as has been found to be the case in the thermodynamic 
limit 0, H3|- This suggests that in general, specific heat 
is a more difficult quantity to calculate numerically. 




0.0 0.5 1.0 1.5 



T/t 

FIG. 12: (Color online) Specific heat as a function of the 
temperature for different densities (from top to bottom n = 
0.7,0.75, . . . , 1.0) and J/t = 0.3. Our Wynn extrapolations 
after four cycles of improvement are compared with FTLM 
results by JP @|. 



VI. CONCLUSIONS 

We have presented an application of the recently in- 
troduced NLC approach to itinerant models. In partic- 
ular, we have studied thermodynamic properties of the 
t-J model on the square lattice. 

From the possible NLC expansions discussed in Ref. 
Q, we have found that the one best suited to the t-J 
model on the square lattice is the site based expansion. 
For this expansion we have shown that the NLC bare 
sums (for the observables considered here) converge to 
lower temperatures than the bare HTE sums. In addi- 
tion, Wynn extrapolations for the site based expansion 
were found to provide better results at lower tempera- 
tures than the ones obtained with Pade extrapolations 
for HTE. This in spite of the fact that the inversion of 
the dependence (n,T) to (n,T) can be done analytically 
within HTE, while within NLC it is done numerically. 

Through a comparison with results obtained by Jaklic 
and Prelovsek using the finite-temperature Lanczos 
method [1, [TTJ] , we have also shown that NLC allows 
one to access regions at low temperatures where finite 
size effects are relevant to exact diagonalization studies. 
One particularly striking example presented here is the 
peak in the specific heat of the Heisenberg model. In the 
calculations done in Refs. [1, [TTJ] it was found that with 
increasing the system size (within the sizes that can be 
addressed by exact diagonalization) the peak does not 
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move while it does increase in size. NLC results, which 
confirmed previous results by Bernu and Misguich [20j . 
exhibit a smaller peak shifted towards lower tempera- 
tures. 

Moving away from Heisenberg models, it remains to be 
seen as to which method is in general accurate to lower 
temperatures. Future application of NLC to supercon- 
ducting and other exotic susceptibilities should prove in- 
formative. The study of the former using HTE remains 
controversial [2l|, [2^, L23[ . It may also prove useful to com- 
bine NLC with Lanczos and FTLM methods to extend 



it to still lower temperatures. 
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